//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// g++ -I /home/tiago/astro_local/lib/ -O2 cur_luz.cpp -o cur_luz
//
// cur_luz.cpp - Nova versao do programa curluz2.c agora escrito em C++.
// Inclui o biblioteca de leitura de arquivos dataFileReader.h (Zabot A.)
// de modo que comentarios podem ser inseridos nos arquivos de entrada. Os
// comentarios devem ser precedidos de #.
// 
// Modo de utilizacao do programa:
// cur_luz arq_in arq_out
// onde arq_in o arquivo com os parametros de entrada e arq_out o arquivo 
// para serem escritos os resultados.
//
// Formato do arquivo de entrada:
// 1a. linha de dados (sem contar linha de comentario):
//   i  - inclinacao (graus).
//   q  - razao de massa (M2/M1).
//   f2 - fluxo da secundaria em fase 0. 
//   fd - fluxo do disco (constante somada a curva de luz).
//   S/N- Relacao Sinal-Ruido na curva
// 2a. linha de dados:
//   limb_law  - lei de obscurecimento de borda.
//   grav_coef - coeficiente do obscurecimento gravitacional.
//   Tpole     - Temperatura do polo da estrela.
//   nspot     - Numero de spots 
// 3a. linha de dados:
//   (coeficientes de obscurecimento de borda)
// 4a. a (4+nspot) linhas de dados: 
//   theta - posicao theta do spot
//   phi   - posicao phi do spot
//   dt    - largura a meia altura em theta
//   dp    - largura a meia altura em phi
//   Is    - Intensidade em relacao a I_L1
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#include <iostream>
#include <cmath>
#include <fstream>
#include <cmath>
#include <vector>
#include <sstream>
#include <string>
#include <cstring>
#include "light_curve.h"
#include "dataFileReader_v0.1.1/dataFileReader.h"

//#include "light_curve.cpp"
//#include "roche.cpp"
//#include "lib_zabot/dataFileReader.h"

using namespace std;

int definicoes_obscurecimento(const int law);
void show_par(ostream& out, vector<vector<double> > &par);

int main(int argc, char *argv[])
{
  // Testar argumentos de entrada...

  if(argc < 3){
    cout << "Modo de utilizacao do programa:\n\n"
	 << "\timOGL arq_p arq_out arq_img\n\n"
	 << "Programa que passa imagens .dat para .ogl\n"
      	 << "Onde arq_p o arquivo com os parametros de entrada e arq_out o\n" 
	 << "arquivo para ser escrita a imagem final. Opcionalmente pode\n"
	 << "se incluir um arquivo com as fases desejadas para ser\n"
	 << "calculados o modelo ([fase.dat]).\n\n"
	 << "Formato do arquivo de entrada:\n"
	 << "1a. linha de dados (sem contar linha de comentario):\n"
	 << "  i  - inclinacao (graus).\n"
	 << "  M1 - Massa da Primaria.\n"
	 << "  M2 - Massa da Secundaria. (q = M2/M1)\n"
	 << "  P  - Periodo Orbital (hrs).\n"
	 << "2a. linha de dados:\n"
	 << "  f2 - fluxo da secundaria em fase 0.\n"
	 << "  grav_coef - coeficiente do obscurecimento gravitacional.\n"
	 << "  limb_law  - lei de obscurecimento de borda.\n"
	 << "  nspot     - Numero de spots.\n"
	 << "  Tpole     - Temperatura do polo da estrela.\n"
	 << "  fd - fluxo do disco (constante somada a curva de luz).\n"
	 << "3a. linha de dados:\n"
	 << "  (coeficientes de obscurecimento de borda)\n"
	 << "4a. a (4+nspot) linhas de dados:\n"
	 << "  theta - posicao theta do spot\n"
	 << "  phi   - posicao phi do spot\n"
	 << "  dt    - largura a meia altura em theta\n"
	 << "  dp    - largura a meia altura em phi\n"
	 << "  Is    - Intensidade em relacao a I_L1\n\n"
	 << "Tiago Ribeiro - OpenGL/LinuxC++ - 2007\n";
    exit(1);
  }

  // Input File Definition - IFD:  
  //     Number of fixed Parameters lines      -  IFD_NP
  //     Number of columns Parameters          - *IFD_P (array[IFD_NP])
  //     Spot Line(onde com. def. de spot - 1) -  IFD_SL
  //     Number of Spot Parameters             -  IFD_SP
  //     posicao da lei de obs. de Borda       -  IFD_B

  const int IFD_NP = 2;
  const int IFD_P[IFD_NP] = {4, 7};
  const int IFD_SL = 3;
  const int IFD_SP = 5;
  const int IFD_B = 2;

  DataFileReader fin(argv[1]);

  vector<vector<double> > in_par(IFD_NP+1);

  //Lendo parametros orbitais

  for(int i = 0; i < IFD_NP; i++){
    ++fin;
    if(fin.numberOfColumns() == IFD_P[i]){
      double data_in;
      for(int j = 0; j < IFD_P[i]; j++){
	data_in = fin.readColumn<double>(j+1);
	if(i == 0 && j == 0) data_in *= M_PI/180;
	in_par[i].push_back(data_in);
      }
    }
    else { 
      ostringstream erro;
      erro << "Numero de parametros errado\n\tArquivo: " 
	   << argv[1] << "\n\tLinha  : " << fin.lineNumber();
      aborta(erro.str());
    }
    /*if(fin.eof()) {
      ostringstream erro;
      erro << "Numero de parametros errado\n\tArquivo: " 
	   << argv[1] << "\n\tLinha  : " << fin.lineNumber();
      aborta(erro.str());
    }*/
  }

  //Lendo coeficientes do obscurecimento de borda
  ++fin;
  if(fin.numberOfColumns()==definicoes_obscurecimento(int(in_par[1][IFD_B]))){
    double data_in;
    for(int i = 0; i < definicoes_obscurecimento(int(in_par[1][IFD_B])); i++){
      data_in = fin.readColumn<double>(i+1);
      in_par[2].push_back(data_in);
    }
  }
  else {
    ostringstream erro;
    erro << "Numero de parametros errado\n\tArquivo: " 
	 << argv[1] << "\n\tLinha  : " << fin.lineNumber();
    aborta(erro.str());
  }


 /* if(fin.eof()) {
    ostringstream erro;
    erro << "Numero de parametros errado\n\tArquivo: " 
	 << argv[1] << "\n\tLinha  : " << fin.lineNumber();
    aborta(erro.str());
  }
  */
  //Lendo parametros do spot

  for(int i = 0; i < int(in_par[1][3]); i++){
    ++fin;
    if(fin.numberOfColumns() == IFD_SP ){
      double data_in;
      in_par.push_back(vector<double> ());
      for(int j = 0; j < IFD_SP; j++){
	data_in = fin.readColumn<double>(j+1) ;
	in_par[IFD_SL+i].push_back(data_in);
      }
    }
    else {
      ostringstream erro;
      erro << "Numero de parametros errado\n\tArquivo: " 
	   << argv[1] << "\n\tLinha  : " << fin.lineNumber();
      aborta(erro.str());
    }

  }

  fin.close();

  show_par(cout, in_par);

  surface SecIm;

  if(argc > 3){
    DataFileReader fin2(argv[3]);

    while(!++fin2){
      if(fin2.numberOfColumns() != 12){
        ostringstream erro;
        erro << "Numero de colunas no arquivo imagem " << argv[3] 
	     << " diferente do esperado!\n"
	     << "\tLido: " << fin2.numberOfColumns() << endl
	     << "\tEsperado: 12" << endl;
        aborta(erro.str());
      }
      double din;
      din = fin2.readColumn<double>(1);
      SecIm.theta.push_back(din);
      din = fin2.readColumn<double>(2);
      SecIm.phi.push_back(din);
      din = fin2.readColumn<double>(3);
      SecIm.px.push_back(din);
      din = fin2.readColumn<double>(4);
      SecIm.py.push_back(din);
      din = fin2.readColumn<double>(5);
      SecIm.pz.push_back(din);
      din = fin2.readColumn<double>(6);
      SecIm.nx.push_back(din);
      din = fin2.readColumn<double>(7);
      SecIm.ny.push_back(din);
      din = fin2.readColumn<double>(8);
      SecIm.nz.push_back(din);
      din = fin2.readColumn<double>(9);
      SecIm.I.push_back(din);
      din = fin2.readColumn<double>(10);
      SecIm.vx.push_back(din);
      din = fin2.readColumn<double>(11);
      SecIm.vy.push_back(din);
      din = fin2.readColumn<double>(12);
      SecIm.vz.push_back(din);
      //fin2 >> din;
    //++fin2;
    
    }
  }

  Roche Star2(SecIm, in_par);
  if(argc == 3) Star2.faz_super();
  //Star2.spot_maker();

  surface SecIm2;
  int size = 100;
  Star2.openGL_Image(size, SecIm2);

  ofstream saida(argv[2]);

  if(!saida){
    ostringstream erro;
    erro << "Arquivo " << argv[2] << " nao pode ser aberto!";
    aborta(erro.str());
  }

  show_par(saida, in_par);
  
  int Flag_pix = 0;
 
  for(int i = 0; i < SecIm2.theta.size(); i+=2){
  //for (int i = 0; i < size; i++){
  // for (int j = 0; j < size*2;j+=2){
      saida << SecIm2.I[i] << " ";
      
    if(Flag_pix >= size) {
      Flag_pix = -1;
      saida << endl;
      //cout << "-> ";
    }
    Flag_pix++;
    //cout << i << endl; //<<i*2*size+j << " " << i << " " << j <<endl;
  }

  saida.close();
  

  return 0;
}

int definicoes_obscurecimento(const int law)
{
  vector<vector<double> > par;

  Light_Curve teste(par);
  
  return teste.npar_obs_borda(law);
}

void show_par(ostream& out, vector<vector<double> > &par)
{
  vector<vector<string> > par_def (4);

  par_def[0].push_back("i = ");
  par_def[0].push_back("M1 = ");
  par_def[0].push_back("M2 = ");
  par_def[0].push_back("P = ");

  par_def[1].push_back("f2 = ");
  par_def[1].push_back("grav_coef = ");
  par_def[1].push_back("limb_law = ");
  par_def[1].push_back("nspot = ");
  par_def[1].push_back("T_pole = ");
  par_def[1].push_back("ffac = ");
  par_def[1].push_back("fd = ");

  for(int i = 0; i < par[2].size(); i++){
    ostringstream str0;
    str0 << "x" << i << " = ";
    par_def[2].push_back(str0.str());
  }

  par_def[3].push_back("s_theta = ");
  par_def[3].push_back("s_phi  = ");
  par_def[3].push_back("dst = ");
  par_def[3].push_back("dsp = ");
  par_def[3].push_back("Is = ");

  out << "# Tiago Ribeiro - LinuxC++ - 2006\n#";

  for(int i = 0; i < par.size(); i++){
    for(int j = 0; j < par[i].size(); j++){
      if( i < 4 ) out << par_def[i][j];
      else out << par_def[3][j];
      
      if( i == 0 && j == 0) out << par[i][j]*180/M_PI << ' ';
      else out << par[i][j] << ' ';
    }
    out << "\n#";
  }
  out << "fase  fluxo  err\n";
}
